#EVALUATION OF DIAGNOSTIC ACCURACY OF MUST & PGSGA

rm(list=ls())

setwd("C:\\Users\\Carol\\Desktop\\Analysis")

if(!R.Version()$arch=="i386"){message("OpenBUGS runs on 32 bit R"); q()}

if(!"BRugs" %in% installed.packages()){install.packages("BRugs")}; require("BRugs")

###################################################################################################

#Conditionally independent Bayesian model - MUST & PGSGA

model.cid <- function(){
  
  #Priors - Test 1: MUST; Test 2:PGSGA
  for (i in 1:2) {
    
    se[i] ~ dbeta(1,1)
    sp[i] ~ dbeta(1,1)
    
  }
  
  for (i in 1:2){
    
    p[i] ~ dbeta(1,1)
    
    pop[i,1:4] ~ dmulti(par[i,1:4],n[i])
    par[i,1] <- se[1]*se[2]*p[i] + (1-sp[1])*(1-sp[2])*(1-p[i])
    par[i,2] <- se[1]*(1-se[2])*p[i] + (1-sp[1])*sp[2]*(1-p[i])
    par[i,3] <- (1-se[1])*se[2]*p[i] + sp[1]*(1-sp[2])*(1-p[i])
    par[i,4] <- (1-se[1])*(1-se[2])*p[i] + sp[1]*sp[2]*(1-p[i])
    
    n[i] <- sum(pop[i,1:4])
    
  }
  
  for(i in 1:2){ #gives PPV & NPV for test 1 and 2 respectively in population 1 (PPV1, NPV1, PPV2, NPV2)
    
    ppv[i] <- p[1]*se[i]/(p[1]*se[i] + (1-p[1])*(1-sp[i]))
    npv[i] <- (1-p[1])*sp[i]/(p[1]*(1-se[i]) + (1-p[1])*sp[i])
    
  }
  
  for(i in 1:2){ #gives PPV & NPV for test 1 and 2 respectively in population 2 (PPV3, NPV3, PPV4, NPV4)
    
    ppv[i+2] <- p[2]*se[i]/(p[2]*se[i] + (1-p[2])*(1-sp[i]))
    npv[i+2] <- (1-p[2])*sp[i]/(p[2]*(1-se[i]) + (1-p[2])*sp[i])
    
  }
  
  for(i in 1:2){
    
    y.index[i] <- (se[i] + sp[i]) - 1 #Youden's index - decide which is best cut.off
    
  }
  
}

#####################################################################################################
#Data

mal.dat <- read.csv("must_pgsga.csv",header=T); mal.dat[,(3:9)] <- lapply(mal.dat[,(3:9)],factor)

str(mal.dat)

cut.offs <- list(c(1,2),c(1,4),c(1,9)); all.results <- list()

for(i in 1:length(cut.offs)){
  
  message(paste('Running model for cut.off',i,sep=" "))
  
  mal.mat <- matrix(NA,ncol=4,nrow=length(unique(mal.dat$site)))
  
  for (p in 1:length(unique(mal.dat$site))){
    
    mal.mat[p,1] <- nrow(mal.dat[mal.dat$site==p & mal.dat$MUST>=cut.offs[[i]][1] & mal.dat$PG.SGA>=cut.offs[[i]][2],]) 
    mal.mat[p,2] <- nrow(mal.dat[mal.dat$site==p & mal.dat$MUST>=cut.offs[[i]][1] & mal.dat$PG.SGA<cut.offs[[i]][2],]) 
    mal.mat[p,3] <- nrow(mal.dat[mal.dat$site==p & mal.dat$MUST<cut.offs[[i]][1] & mal.dat$PG.SGA>=cut.offs[[i]][2],]) 
    mal.mat[p,4] <- nrow(mal.dat[mal.dat$site==p & mal.dat$MUST<cut.offs[[i]][1] & mal.dat$PG.SGA<cut.offs[[i]][2],]) 
    
  }
  
  mal.list <- list(mal.mat); names(mal.list) <- "pop"
  
  ########################################################################################################
  #Convert all to BUGS format
  
  #write model to a file
  writeModel(model.cid,'model.cid.txt')
  
  #Bugs data
  bugsData(mal.list,fileName='mal.list.txt')
  
  #make 2 initial values chains 
  bugsInits(inits=list(list(se=rep(0.85,times=2),sp=rep(0.60,times=2),p=rep(0.05,times=2))),numChains=1,'CID.Init1.txt')
  bugsInits(inits=list(list(se=rep(0.75,times=2),sp=rep(0.85,times=2),p=rep(0.50,times=2))),numChains=1,'CID.Init2.txt')
  
  #now check, load data, compile etc.
  modelCheck("model.cid.txt") #check model file.
  modelData("mal.list.txt") #read data file
  modelCompile(numChains=2) #compile model with 2 chains
  modelInits('CID.Init1.txt',1) #read init data file
  modelInits('CID.Init2.txt',2) #read init data file
  
  #modelGenInits() #generate the missing initial values
  
  modelUpdate(3500) #burn in
  
  samplesSet(c('se','sp','p','ppv','npv','y.index')) #parameters to monitor
  
  modelUpdate(6000) #more iterations 
  
  #SOME DIAGNOSTICS FIRST
  #Check convergence (Trace plots) - should ideally check all
  #samplesHistory('se',mfrow=c(1,1)) # plot the chain,
  #samplesHistory('sp',mfrow=c(1,1)) # plot the chain,
  #samplesHistory('p',mfrow=c(1,1))
  
  #Plot the Gelman-Rubin diagnostic statistics - ratio should be close to 1
  #samplesBgr('se',mfrow=c(1,1))  
  #samplesBgr('sp',mfrow=c(1,1))
  #samplesBgr('p',mfrow=c(1,1))
  
  #Density plots
  #samplesDensity('se',mfrow=c(1,1)) 
  #samplesDensity('sp',mfrow=c(1,1))
  #samplesDensity('p',mfrow=c(1,1))
  
  (all.results[[i]] <- samplesStats('*')); names(all.results)[i] <- paste('cut.off_',i,sep="")
  
}

capture.output(all.results,file="all.results.txt")

##################################################################################################
